Improved Method for Accurate and Eecient Quantiication of Mrs Data with Use of Prior Knowledge

We introduce AMARES (Advanced Method for Accurate, Robust and E cient Spectral tting), an improved method to accurately and e ciently estimate the parameters of noisy magnetic resonance spectroscopy (MRS) signals in the time domain. As a reference time domain method we take VARPRO. VARPRO uses a simple Levenberg-Marquardt algorithm to minimise the variable projection functional. This variable projection functional is derived from a general functional, which minimises the sum of squared di erences between the data and the model function. AMARES minimises the general functional which improves the robustness of MRS data quanti cation. The newly developed method uses a version of NL2SOL, a sophisticated nonlinear least squares algorithm, to minimise the general functional. In addition, AMARES uses a singlet approach for imposition of prior knowledge instead of the multiplet approach of VARPRO because this greatly extends the possibilities of the kind of prior knowledge that can be invoked. Other new features of AMARES are the possibility of tting echo signals, choosing a Lorentzian as well as a Gaussian line shape for each peak and imposing lower and upper bounds on the parameters. Simulations, as well as in-vivo experiments, con rm the better performance of AMARES compared to VARPRO in terms of accuracy, robustness and exibility. 2 Introduction For medical diagnosis or biochemical analysis accurate and e cient quanti cation of magnetic resonance spectroscopy (MRS) signals is of utmost importance. MRS signals however are often characterised by a low signal-to-noise ratio and overlapping peaks. Under these circumstances simple signal processing algorithms like numerical integration are not adequate. In recent years a number of more advanced time domain techniques based on a mathematical model function have been developed. Noninteractive methods exist that are noniterative, computationally e cient and which can be fully automatic. A serious drawback however is the fact that only very limited prior knowledge can be incorporated in these algorithms and that the model function is restricted to a sum of complex exponentially damped sinusoids. Among this class of methods are the algorithms based on Kumaresan-Tufts' linear prediction (LP) method (1) combined with the singular value decomposition (SVD) (2). Kung's state-space approach (3) combined with SVD (called HSVD (4)) is a more e cient and more accurate alternative to the LP methods as it circumvents the polynomial rooting and root selection. Rapid and more accurate variants of the state-space algorithms have been recently proposed (5{8), but the limitations to the imposition of prior knowledge about model function parameters are inherent to this type of methods. On the other hand, interactive methods exist that are iterative, with more user-involvement, in some cases computationally less e cient, but that do allow inclusion of prior knowledge. The algorithms t the data to the nonlinear model function in a least squares sense, leading to maximum likelihood parameter estimates in the case of white Gaussian noise. See (9) for an overview of time domain methods. In this paper we introduce AMARES, an improved interactive time domain method for accurate and e cient parameter estimation of MRS signals with use of prior knowledge. As a reference interactive method we take VARPRO (10), a method that has proven to be very reliable 3 in recent years (11,12). Our method improves VARPRO in several ways. First of all, we minimise a di erent functional than in VARPRO leading to more robust quanti cation of MRS signals. In addition, we use a more sophisticated nonlinear least squares (NLLS) algorithm which has a more robust behaviour than the simple Levenberg-Marquardt (LM) algorithm used in VARPRO. The robustness is also increased by imposing lower and upper bounds on the parameters. This allows to impose the physical requirement that all damping factors are positive. It is shown here that this strategy leads to greater robustness than the transformation of variables used in VARPRO. We also greatly extended the possibilities of imposing prior knowledge. This leads to more accurate parameter estimation since the precision of the quantitative analysis is increased by using all available information (9) and increases the exibility for the user. In addition, using all available prior knowledge leads to an optimisation problem in a smaller number of variables which can be solved more consistently in less time. As extra features of our new algorithm we made it possible to use a more general model function than in VARPRO and to t entire spin echo signals instead of having to truncate the rst part of the echo as in VARPRO. As a result our method allows for more robust, accurate, exible and e cient quanti cation of MRS signals. In the rst part of the paper we will explain the algorithmic di erences between VARPRO and AMARES and we will point out the extended possibilities of imposing prior knowledge in AMARES. The di erences in robustness and e ciency between VARPRO and AMARES are illustrated in the second part of the paper. In the third part we illustrate the extended possibilities of imposing prior knowledge using an in vivo signal. 4 Methods The most commonly used model function in MRS signal processing to t N measured data points yn is a sum of exponentially damped sinusoids (Lorentzian lines after FT): yn = ŷn + en = K Xk=1 akej ke( dk+j2 fk)tn + en n = 0; 1; : : : ; N 1 [1] where j = p 1, ak is the amplitude, k the phase, dk the damping factor and fk the frequency of the kth sinusoid (k = 1; : : : ;K); tn = n t+ t0 with t the sampling interval (non-uniform sampling vectors are also valid), t0 the time between the e ective time origin and the rst data point to be included in the analysis and en is complex white Gaussian noise. The caret on y indicates that this quantity represents the model function rather than the actual measurements. Other types of model line forms can be used; we will come back to this. To obtain maximum likelihood estimates in the case of white Gaussian noise, one has the choice between minimising two functionals. The rst one, denoted by G throughout the paper, is straightforwardly derived using probability theory and is given by the following formula: G(a;d; f ; ; t0) = N 1 X n=0 jyn K Xk=1 akej ke( dk+j2 fk)tn j2 = ky lk2; [2] with y=[y0; : : : ; yN 1]T the signal vector, l =[a1ej 1 ; : : : ; aKej K ]T , a;d; f ; are the vectors of amplitudes, dampings, frequencies and phases respectively. The superscript T denotes the 5 transpose, k:k the Euclidean vector norm and = 26666664 e( d1+j2 f1)t0 : : : e( dK+j2 fK)t0 ... . . . ... e( d1+j2 f1)tN 1 : : : e( dK+j2 fK)tN 1 37777775 [3] is an N K matrix of full rank. The second functional is derived as follows. Suppose the nonlinear parameters f and d are known, then the matrix can be computed and an estimate for the linear parameters l is obtained by solving a linear LS problem: bl = yy, with y = ( T ) 1 T , the pseudo-inverse of . Substituting this estimate into the original functional G of Eq. [2] results in: V (d; f ; t0) = ky yyk2 [4] which is called the variable projection functional (13) and denoted by V from now on. In this way the amplitudes a and the phases are eliminated. As a result we obtain a minimisation problem where the number of variables is reduced but where the functional has become more complicated. Both functionals consist of a sum of squared residuals and give rise to typical NLLS problems. The available biochemical prior knowledge can be expressed as a set of linear relations between parameters resulting in a minimisation problem with linear equality constraints. These constraints are substituted in the original functional in order to obtain an unconstrained NLLS minimisation problem. As a consequence, regardless which of the two functionals in Eq. [2] or Eq. [4] is used, we always have to solve an unconstrained NLLS optimisation problem. The latter can be solved using local or global optimisation theory. Global optimisation has already been used in MRS (14, 15) but the main disadvantage of these methods is the poor computational e ciency. The methods are too time consuming to be useful in practical 6 situations. On the other hand it is possible to obtain good starting values by peak picking for use in local optimisation methods, resulting in an acceptable solution in a reasonable time in most circumstances. Therefore, in this paper we focus our attention on methods based on local optimisation theory. VARPRO: algorithmic aspects VARPRO minimises the variable projection functional V using a modi ed version of Osborne's Levenberg-Marquardt algorithm (16). The basic model function of Eq. [1] is extended so that one has the choice between a pure Gaussian (g=1) or a pure Lorentzian (g=0) line shape for the entire signal: yn = ŷn + en = K Xk=1 akej ke( b2k(1 g+gtn)tn)ej2 fktn + en n = 0; 1; : : : ; N 1 [5] It is very important to note here that in VARPRO minimisation is done w.r.t. the square roots of the damping factors bk, i.e., b2k = dk (k = 1; : : : ;K), to make sure the damping is positive as required by physical considerations. This variable transformation works in most practical situations, but is not without danger as explained by Gill and co-workers (17) and as we will show in the section \Simulations". Di erent forms of multiplet prior knowledge can be imposed as relations between parameters of the same type, a detailed explanation is given in the subsection \VARPRO: multiplet approach for prior knowledge". With VARPRO it is possible to do frequency-selective quanti cation in the time domain (18). As a measure of precision we compute approximations to the Cram er-Rao lower bounds as explained in (9). 7 AMARES: algorithmic aspects AMARES uses dn2gb (available from the Port library of netlib (19)), the most recent version of NL2SOL to minimise the general functional. NL2SOL is a secant method recommended in (20) and has the advantage of handling large-residual or very nonlinear problems better than a Levenberg-Marquardt algorithm (as used in VARPRO). An additional advantage of using dn2gb is the fact that the algorithm allows the user to specify upper and lower bounds on the variables. We can use this feature to impose the natural bounds on the variables. The physical requirement of a positive damping can thus easily be imposed. Likewise amplitudes (like the dampings) cannot become negative (negative amplitudes correspond to a 180 phase shift), whereas the upper and lower bounds on the frequencies are determined by the spectral width. Phases can be constrained to lie between 180 and +180 . In general it is recommended to add those extra bounds on the variables in order to ensure maximum accuracy and robustness (17). We show in the section \Simulations" that this new approach of ensuring the positiveness of the dampings leads to a more robust method than VARPRO. AMARES ts the data to the following model function in a least squares sense: yn = ŷn + en = K Xk=1 akej ke( dk(1 gk+gktn)tn)ej2 fktn + en n = 0; 1; : : : ; N 1 [6] The basic model function of Eq. [1] is extended so that one now has the choice between a pure Lorentzian (gk = 0) or a pure Gaussian line form (gk = 1) for every peak separately as opposed to VARPRO where the line form used has to be the same for all peaks. AMARES allows the imposition of all kinds of linear relations between individual parameters by using a singlet instead of a multiplet approach, thereby greatly increasing the prior knowledge that can be used in the data processing. This is explained in more detail in the 8 subsection \AMARES: singlet approach for prior knowledge" and is illustrated in the section \Quanti cation of an in vivo signal using prior knowledge". AMARES also o ers the possibility to t echo signals, an echo being modeled as two FIDs back to back. The left and the right part of the echo are considered to have the same amplitudes, frequencies and phases but di erent dampings. The dampings of the right and left part can however be linked to each other. In VARPRO the rst part of the echo is truncated in order to be able to work with the model function of a FID signal. As a consequence, part of the signal and thus part of the information is destroyed. Like in VAPRRO the user has the possibility of doing frequency-selective quanti cation in the time domain (18). The approximations to the Cram er-Rao lower bounds are used as a measure of precision. Note that the Cram er-Rao bounds do not change by imposing simple bounds on the variables as proven in (21). Computational considerations VARPRO and AMARES both make use of evaluations of the residuals and the Jacobian, the Jacobian being the matrix consisting of the rst derivatives of the residuals with respect to (w.r.t.) the unknown parameters. In both methods the residuals and Jacobian are provided analytically. Evaluation of the residuals and the Jacobian of V (VARPRO) is far more complicated than in the case of G (AMARES). In VARPRO a simpli cation introduced by Kaufman (22) is used which results in a more e cient computation of the Jacobian of V . To perform the actual minimisation, the data and the model function are split in a real and imaginary part as follows. V = kyV V lV k = kyV V Vyyk [7] 9 where V = 666666664 e d1t0cos(2 f1t0) e d1t0sin(2 f1t0) : : : e dKt0cos(2 fK t0) e dKt0sin(2 fKt0) e d1t0sin(2 f1t0) e d1t0cos(2 f1t0) : : : e dKt0sin(2 fK t0) e dKt0cos(2 fK t0) ... ... . . . ... ... e d1tN 1cos(2 f1tN 1) e d1tN 1sin(2 f1tN 1) : : : e dKtN 1cos(2 fK tN 1) e dKtN 1sin(2 fKtN 1) e d1tN 1sin(2 f1tN 1) e d1tN 1cos(2 f1tN 1) : : : e dKtN 1sin(2 fK tN 1) e dKtN 1cos(2 fK tN 1) 777777775 [8] yV = 66666666666664 Re(y0) Im(y0) ... Re(yN 1) Im(yN 1) 77777777777775 lV = 66666666666664 a1cos( 1) a1sin( 1) ... aKcos( K) aKsin( K) 77777777777775 = 66666666666664 Re(c1) Im(c1) ... Re(cK) Im(cK) 77777777777775 [9] with ck = akej k ; (k = 1; : : : K). G = kyG GlGk [10] where G = 66666666666664 e d1t0cos(2 f1t0 + 1) : : : e dK t0cos(2 fKt0 + K) e d1t0sin(2 f1t0 + 1) : : : e dK t0sin(2 fKt0 + K) ... . . . ... e d1tN 1cos(2 f1tN 1 + 1) : : : e dK tN 1cos(2 fKtN 1 + K) e d1tN 1sin(2 f1tN 1 + 1) : : : e dK tN 1sin(2 fKtN 1 + K) 77777777777775 [11] 10 yG = 66666666666664 Re(y0) Im(y0) ... Re(yN 1) Im(yN 1) 77777777777775 lG = 26666664 a1 ... aK 37777775 [12] Re(:), Im(:) denote the real, respectively imaginary part of a complex quantity. In VARPRO and AMARES starting values for the frequencies and dampings need to be provided by the user. In AMARES, the starting values for the amplitudes and phases are computed by solving the LS problem blV = V yyV , with the starting values for the frequencies and dampings inserted in V. VARPRO and AMARES are both written in Fortran77. They can be used as stand-alone programs, but they are also implemented within the MRUI software package (23), a graphical user interface to facilitate use of sophisticated spectral analysis routines in biomedical/biochemical laboratories and the clinic. AMARES is available from the authors upon request. VARPRO: multiplet approach for prior knowledge VARPRO puts peaks belonging to the same multiplet into one group. Amplitudes within one group can be left unconstrained or they can be linked to each other when the intensity ratios of the peaks are known. Dampings are treated in the same way, but with the extra possibility of xing the damping to a value which has to be the same for all peaks in that group. Frequencies within a group can be left unconstrained or they can be linked in case the frequency splittings between neighbouring peaks are known exactly. If the frequency splittings between neighbouring peaks are unknown, but equal for all the peaks in the group, a new variable is introduced. One then has f1, f2 = f1+ , f3 = f1+2 , : : : , fk = f1+(k 1) , for all k peaks belonging to the same group, fk is the arbitrarily chosen reference peak, neighbouring peaks in the spectrum 11 are denoted by consecutive numbers. In addition, the frequency of the rst component (chosen arbitrarily) can be xed together with the frequency splitting of the multiplet. In that case all the frequencies of a group are xed. Note that it is impossible to x the frequencies of the peaks of a certain group independently: the frequency di erences between neighbouring peaks are the same. The overall phase k of a peak is considered to consist of an individual phase 0k and a so-called zero order phase 0 which is the same for all peaks: k = 0k + 0. The individual phases can be left unconstrained in which case the zero order phase must be xed. The phases of all peaks within a group can be xed to a value (the same for all the peaks within the group) relative to the zero order phase, which can be estimated or kept xed. Of all the groups with constraints on amplitudes one group can be taken as a reference and the other groups can be linked to the reference one. This means that the relative amplitude of the rst component (chosen arbitrarily) of the reference group w.r.t. the rst component of a linked group can be imposed. Groups with constrained dampings are treated in the same way. Groups with variable frequency splittings can also be linked in the same way with as only di erence that the frequency splitting ratio between the linked groups has to be one. We illustrate some of the above mentioned possibilities using the 31P simulation example. The peaks belonging to the -ATP triplet are put into one group. The amplitude ratios are a1 : a2 : a3 = 1 : 2 : 1. If we impose this prior knowledge in VARPRO and we leave the corresponding phases unconstrained, the program will automatically assume these phases to be equal (this is not the case when amplitudes are left unconstrained). This is equivalent to linking the complex amplitudes akej k or in this example Re(c1) : Re(c2) : Re(c3) = 1 : 2 : 1 and Im(c1) : Im(c2) : Im(c3) = 1 : 2 : 1. Peaks 4 and 5 belong to the -ATP doublet and are put in a second group. We know that a4 : a5 = 1 : 1 and that a4 : a1 = 2 : 1. VARPRO o ers the possibility of imposing this prior knowledge between the two groups as long as the constraints on the phases are equal for both groups. 12 In the case of the -ATP triplet we can impose a xed frequency splitting of 16 Hz between the individual peaks within the triplet, leading to: f2 = f1 + 16 and f3 = f1 + 32. Suppose however that we do not know the exact frequency splittings between the peaks but that we do know that these splittings are the same. We can express this by introducing a new variable . We then get: f2 = f1 + and f3 = f1 + 2 . We can also express that the frequency splitting in the -ATP doublet is the same as in the -ATP triplet. Using VARPRO we express this as f5 = f4 + . This implementation of prior knowledge in VARPRO already o ers a lot of possibilities of imposing various linear relations between parameters and has proven crucial in earlier studies (12,24,25). However, in recent years, more prior knowledge has become available and the present implementation of VARPRO can no longer satisfy all the needs (11). One example is formed by the six glycogen resonances in 13C MRS, for which the relative frequency shifts are known but di erent. This can therefore not be implemented as a regular multiplet in the VARPRO algorithm. For the quanti cation of a 13C MRS signal using AMARES, we refer to (26). AMARES: singlet approach for prior knowledge In AMARES we have chosen for a singlet approach making an identical treatment of all parameters possible. This approach allows everything that was previously possible using the multiplet approach of VARPRO and much more. Each of the parameters can be left unconstrained or kept xed. One can impose a xed shift or ratio w.r.t. any unconstrained parameter of the same type. Finally, it is possible to impose a variable shift or ratio w.r.t. any unconstrained or xed parameter of the same type. These variable shifts or ratios can then be linked between di erent groups of peaks. Imposing a xed shift or ratio and a variable shift or ratio deserves some more explanation, the other forms of prior knowledge are trivial to implement. As an arbitrary example, take a signal consisting of three peaks. Imposing a xed ratio is 13 done in the following way; e.g. for the amplitudes: a2 = a1x; a3 = a1y; x; y 2 R. x and y are values provided by the user. Imposing a variable ratio leads to the introduction of a new variable , a2 = a1(x ); a3 = a1(y ); x; y 2 R. x and y are values provided by the user. Imposing a variable shift or ratio is only useful if the new variable is shared by at least two peaks, otherwise the total number of unknowns is not changed. Simulations In this section we demonstrate the di erences between VARPRO and AMARES. To this end we perform a Monte Carlo study. The simulation signal we used was derived from an in vivo 31P spectrum measured in the human brain and consisted of 256 complex data points and 11 exponentials (see Figure 1 and Table 1). The 31P peaks from brain tissue, from left to right phosphomonoesters, inorganic phosphate, phosphodiesters, phosphocreatine, -ATP, -ATP and -ATP can all be observed. From the noiseless signal 300 noisy realisations were generated with noise standard deviation (both on the real and imaginary parts). We used a low, intermediate and high noise level ( = 5; 15; 25). A method is considered to fail if the minimisation algorithm claims not to have found the solution (according to its corresponding stopping criterium) or if not all 11 peaks are resolved within speci c intervals lying symmetrically around the exact frequencies. The halfwidths of these intervals are 8.6, 7.3, 8.7, 3.2, 3.2, 3.5, 3.6, 0.7, 5.6, 2.4 and 7.8 Hz, respectively. These values are derived from the Cram er-Rao lower bounds on the frequencies at that noise standard deviation where the two intervals of two neighbouring peaks in the -ATP triplet touch. For all our results we veri ed that these simple criteria guarantee that all methods nd the same minimum if no failure occurs. As a result, root mean-squared error, bias and standard deviation of the parameter estimates computed by any method as a function of the noise level are the same and will not be shown. 14 Starting values for the dampings and the frequencies can be obtained by peak picking or by using a noninteractive method like HSVD. In a previous study (27) we found that for signals with a high signal-to-noise ratio HSVD is the best way to provide starting values. However for low signal-to-ratio signals peak picking is preferred. Here only peak picking is considered. For every noise level used in the Monte Carlo simulation we randomly picked one signal out of 300 for which damping and frequency starting values were determined by peak picking. These starting values were then used to process all the signals a ected by the same noise level. We compared VARPRO and AMARES in terms of robustness and e ciency. The robustness could be assessed by the number of times a method fails. To compare the e ciency we looked at the average number of functional and Jacobian evaluations and the average overall CPU time. All experiments were performed on a SUN ULTRA2 (200 MHz). In the rst experiment no prior knowledge is imposed, hence 4 parameters for each of the 11 peaks are to be estimated. AMARES minimises a problem in 44 variables, whereas VARPRO minimises a problem in 22 variables. The results are displayed in Table 2. We see that VARPRO performs poorly in terms of failures compared to AMARES. We did several tests to nd the reason for this signi cantly poorer robustness of VARPRO. First, we changed the VARPRO code so that G was minimised instead of V . The results are displayed in Table 3. We see that the failure rate drops to a level comparable to that obtained with AMARES. This suggests that minimisation of G is easier than of V . We see an equal improvement in failure rate when we take the VARPRO code and modify it so that minimisation is done w.r.t. dk; k = 1; : : : ;K directly. The dampings are forced to be positive by imposing a penalty on the function if the algorithm computes damping estimates that are negative. Note however that this approach can lead to slow convergence in some cases as pointed out in (28). To investigate the in uence of the minimisation algorithm used, we replaced the LM algo15 rithm in VARPRO by the NL2SOL algorithm used in AMARES. The results are displayed in Table 3. Although the functional minimised is V and the minimisation is done w.r.t. bk; k = 1; : : : ;K, the failure rate drops to the same level as the one obtained using AMARES. This clearly shows that AMARES is a more robust algorithm, able to deal with more di cult optimisation problems. When more prior knowledge is imposed or the starting values improve, leading to a reduction in complexity of the optimisation problem, the di erences in robustness between VARPRO and AMARES decrease. In these situations VARPRO can still be more e cient, due to the less complicated nature of the underlying optimisation algorithm. Quanti cation of an in vivo signal using prior knowledge To illustrate the features of the new algorithm we tested it on a challenging data set provided by Dr. A. Heerschap of the University Hospital Nijmegen (see acknowledgments). An in vivo spectrum was acquired from the prostatic gland of a benign prostatic hyperplasia patient, at 1.5T using a PRESS localisation sequence ( 1=11ms, 2=67.5 ms) and an endorectal coil for signal reception. 256 scans of 1024 complex FID data points were acquired from a volume of 8 cc, with a TR of 1.6 s (see Figure 2). The goal of this measurement was the quanti cation of the citrate content, which is an important tool in the discrimination between prostate adenocarcinoma and benign prostatic hyperplasia (BPH). The citrate protons form a complicated AB-type multiplet, consisting of four Lorentzian lines with di erent phase. Relative intensities of the citrate lines and their phase deviations from the rest of the spectrum can be calculated (29) and were provided to us by Dr. M. van der Graaf (see acknowledgments). It is also known that the linewidths of the main peaks should be equal as well as those of the side peaks. To analyse the signal with VARPRO we put all the peaks in separate groups. We could 16 have put the peaks of choline and creatine in one group but that would not have made any di erence. The four peaks of citrate have to be put in separate groups in order to impose the xed (but di erent for each peak) phase shifts w.r.t. the zero order phase. This implies that we cannot impose the known frequency shifts between the citrate peaks since these can only be imposed between peaks belonging to the same group. Amplitudes were linked between the di erent groups and the individual phases of the peaks xed relative to the zero order phase. Since VARPRO allows only one overall linking for dampings between groups, we could link all four peaks to each other (theoretically incorrect), or either the side peaks, or the main peaks. In this example with the VARPRO method the linewidths of the small side peaks were linked, as these are more di cult to t. Zero order phase and begin time are estimated. The results are displayed in Table 4. Only with the new algorithm we are capable of imposing all the available prior knowledge. The results are displayed in Table 5. The most important di erence between the VARPRO results and those of the new method is that the latter can be obtained in a much more consistent fashion. This is entirely due to the increased use of important prior knowledge on the AB-type multiplet of citrate, in particular the chemical shift di erences between its respective components. This is especially crucial at echo times where the side peaks reach their minimum intensities and sink into the noisy baseline (30, 31). With the combined prior knowledge of the chemical shift di erences and amplitude (peak area) ratios with respect to the prominent main peaks of the quartet, the small side peaks can still be estimated reliably. Including the small side peaks in the t is not only necessary for complete quantitation of the citrate multiplet, but also enhances the accuracy of the creatine t, as their resonances overlap closely (see Figure 3). Another advantage of the new method is that, in combination with all other prior knowledge, the linewidths of the side peaks can be linked to each other, and those of the main peaks can be linked to each other, see the results in Table 5. The main peaks were found to be a little broader than the side 17 peaks, which corresponds to the expectations for the citrate multiplet (32). With the VARPRO method where we could only link the side peaks to each other, we obtained slightly di erent (although not signi cant, due to the high level of the noise) linewidths of the main peaks, which is undesired. Despite the slight di erences in linewidths and chemical shifts between the two methods, the estimated amplitudes are identical for the citrate multiplet. This is mainly due to the available prior knowledge for the amplitude ratios with respect to the prominent main peaks, and to the precise phase constraints. As an aside, it can be noted that, perhaps because of the improved tting of the left side peak of citrate, the tting of creatine and choline is more according to our expectations for the new method. Their linewidths (damping factors) should not di er too much, due to the leveling e ect on the e ective linewidths of the magnetic eld inhomogeneities and magnetic susceptibility e ects. Creatine may indeed be a little broader than choline, as its T2 is shorter than that of choline (33). Conclusions We have presented AMARES, a new algorithm which outperforms the currently used reference time domain method VARPRO in two ways. First of all we made an appropriate choice of the functional and the nonlinear least squares algorithm. When estimating the parameters of a magnetic resonance spectroscopy signal in the time domain using optimisation methods one has the choice between two functionals. On the one hand, we can minimise a general functional, consisting of the sum of squared di erences between the data and the model function. On the other hand, a so-called variable projection functional can be derived from the general functional by eliminating all linear parameters and be optimised. We showed here that minimising the general functional leads to a more robust determination of parameter estimates. In addition we showed that the functional should be minimised with respect to the damping with imposition of a positivity constraint instead of minimising the functional with respect to the square root of 18 the damping as is done in VARPRO. In this way numerical problems are avoided.Secondly, AMARES allows the inclusion of more prior knowledge about the signal param-eters, the model function (Lorentz, Gauss) and the type of signal (FID or echo). This resultsin increased accuracy and userexibility. In addition, imposing all prior knowledge leads to aproblem in a smaller number of variables which can be solved in less time. Overall we obtainan algorithm which outperforms VARPRO in terms of accuracy, robustness and exibility.AcknowledgmentsThe authors thank Dr. Marinette van der Graaf and Dr. Arend Heerschap of the Departmentof Radiology, University Hospital Nijmegen St. Radboud, for supplying the example data setand the relevant molecular prior knowledge.The rst author is a Ph.D. student funded by the IWT (Flemish Institute for Support ofScienti c-Technological Research in Industry). The third author is a Research Associate withthe F.W.O. (Fund for Scienti c Research-Flanders). This work is supported by the Belgian Pro-gramme on Interuniversity Poles of Attraction (IUAP-4/2 & 24), initiated by the Belgian State,Prime Minister's O ce for Science, Technology and Culture, by the EU Programme \HumanCapital and Mobility, Networks", project CHRX-CT94{432, and by a Concerted Research Ac-tion (GOA) project of the Flemish Community, entitled \Model-based Information ProcessingSystems".19 References1 . R. Kumaresan and D.W. Tufts, IEEE Trans. ASSP 30, 833 (1982).2 . H. Barkhuysen, R. de Beer, W.M.M.J. Bovee, and D. van Ormondt, J. Magn. Reson. 61,465 (1985).3 . S.Y. Kung, K.S. Arun, and D.V. Bhaskar Rao, J. Opt. Soc. Am. 73, 1799 (1983).4 . H. Barkhuysen, R. de Beer, and D. van Ormondt, J. Magn. Reson. 73, 553 (1987).5 . W.W.F. Pijnappel, A. van den Boogaart, R. de Beer, and D. van Ormondt, J. Magn. Reson.97, 122 (1992).6 . S. Van Hu el, H. Chen, C. Decanniere, and P. Van Hecke, J. Magn. Reson. A 110, 228(1994).7 . H. Chen, S. Van Hu el, C. Decanniere, and P. Van Hecke, J. Magn. Reson. A 109, 46(1994).8 . H. Chen, S. Van Hu el, D. van Ormondt, and R. de Beer, J. Magn. Reson. A 119, 225(1996).9 . R. de Beer and D. van Ormondt, \NMR Basic Principles and Progress", M. Rudin Ed.,Vol. 26, p. 201, Springer-Verlag, Berlin, 1992.10 . J.W.C. van der Veen, R. de Beer, P.R. Luyten, and D. van Ormondt, Magn. Reson. Med.6, 92 (1988).11 . C. Decanniere, P. Van Hecke, H. Chen, S. Van Hu el, C. van der Voort, B. van Tongeren,and D. van Ormondt, J. Magn. Reson. B 105, 31 (1994).12 . A. van den Boogaart, F.A. Howe, L.M. Rodrigues, M. Stubbs, and J.R. Gri ths, NMRBiomed. 8, 87 (1995).20 13 . G.H. Golub and V. Pereyra, SIAM J. Numer. Anal. 10, 413 (1973).14 . F.S. DiGennaro and D. Cowburn, J. Magn. Reson. 96, 582 (1992).15 . G.J. Metzger, M. Patel, and X. Hu, J. Magn. Reson. B 110, 316 (1996).16 . M.R. Osborne, \Numerical methods for non-linear optimization", Lootsma Ed., AcademicPress, London, 1972.17 . P.E. Gill, W. Murray, and M.H. Wright, \Practical Optimization", Academic Press, SanDiego, 1988.18 . A. Knijn, R. de Beer, and D. van Ormondt, J. Magn. Reson. 97, 444 (1992).19 . ftp:==netlib.att.com.20 . J.E. Dennis and R.B. Schnabel, \Numerical methods for unconstrained optimisation andnonlinear equations", Prentice Hall, Englewood Cli s, 1983.21 . J.D. Gorman and A.O. Hero, IEEE Trans. on Inf. Theory 26, 1285 (1990).22 . L. Kaufman, BIT 15, 49 (1975).23 . A. van den Boogaart, S. Cavassila, L. Vanhamme, J. Totz, and P. Van Hecke, In Proceedings18th Ann. Int. Conf. IEEE Eng. Med. Biol. Soc., Amsterdam, The Netherlands (1996). Seealso http:==mrui-web.uab.es=mruiwww=mrui hom.html.24 . A. van den Boogaart, M. Ala-Korpela, J. Jokisaari, and J.R. Gri ths, Magn. Reson. Med.31, 347 (1994).25 . M. Stubbs, A. van den Boogaart, C.L. Bashford, P.M.C. Miranda, L.M. Rodrigues, F.A.Howe, and J.R. Gri ths, Biochim. Biophys. Acta 1291, 143 (1996).26 . L. Vanhamme, S. Van Hu el and A. van den Boogaart, Proc., ISMRM, pp. 1414, 5thScienti c Meeting, Vancouver, Canada (1997).21 27 . L. Vanhamme, A. van den Boogaart, and S. Van Hu el, Proc. EUSIPCO-96, Trieste, Italy(1996).28 . J.J. More, B.S. Garbow and K.E. Hillstrom, User Guide for Minpack-1, Argonne, IL,ANL-80-74 (1980).29 . M. van der Graaf, A. van den Boogaart,W.M.M.J. Bovee, and A. Heerschap, Proc., ISMRM,pp. 1419, 5th Scienti c Meeting, Vancouver, Canada (1997).30 . M. van der Graaf, G.J. Jager, and A. Heerschap, MAGMA 5, in press (1997).31 . M. van der Graaf, G.J. Jager, and A. Heerschap, Proc. ESMRMB XIII, MAGMA 4 vol. 2suppl., 320 (1996).32 . W.M.M.J. Bovee, J. Magn. Reson. 24, 327 (1976).33 . A. van den Boogaart, \The use of Signal Processing Algorithms to Obtain BiochemicallyRelevant Parameters form Magnetic Resonance Data Sets", Ph.D. thesis, U. London, 1995.22 Table 1: Exact parameters values of the simulated 31P MRS signal, modeled by Eq. [1]. t0 =0; t = 0:333. The values are based on the t of an in vivo 31P MRS signal.peak k fk (Hz) dk(Hz) ak (a.u.)k ( ) a18650751352705015013535450751354152501501355168501501356292501501357308501501358360251501359440285:7140013510490256013511530200500135a k = k 180= expresses the phase in degrees23 Table 2: Sample Mean and standard deviation values of Monte Carlo simulation (300 runs)no prior knowledge. VARPRO and AMARES are compared. The computed quantities are:% fail, the percentage of failures, jev, the average number of Jacobian evaluations, fev, theaverage number of functional evaluations, cpu, the average CPU time in seconds. std denotesthe standard deviation.criterium VARPRO AMARES5% fail00jev std 10:4 1:0 11:5 1:7fev std 13:5 1:3 15:6 2:4cpu std 2:6 0:3 2:1 0:315 % fail12.30.33jev std 11:5 3:8 14:7 3:2fev std 15:8 6:9 20:9 5:8cpu std 2:9 1:0 2:6 0:625 % fail13.07.0jev std 11:3 4:2 15:0 6:1fev std 14:1 5:8 20:4 9:7cpu std 2:8 1:0 2:6 1:124 Table 3: Sample Mean and standard deviation values of Monte Carlo simulation (300 runs)no prior knowledge. VARPRO-G is the original VARPRO code adapted to minimise the gen-eral functional G instead of the variable projection functional V . VARPRO-d is the originalVARPRO code modi ed to minimise w.r.t. dk; k = 1 : : : K. VARPRO-N is the original VARPROcode where the LM algorithm is replaced by NL2SOL. Further as in Table 2.criterium VARPRO-G VARPRO-d VARPRO-N5% fail000jev std 9:9 0:39:0 0:38.5 0.5fev std 10:9 0:312:1 0:511.5 0.5cpu std 1:7 0:12:2 0:072.1 0.115 % fail0.3300jev std 10:6 1:310:2 0:911.8 1.5fev std 11:7 1:714:4 1:114.8 2.0cpu std 1:8 0:22:8 0:22.7 0.325 % fail7.06.76.7jev std 12:7 4:511:4 3:512.3 3.7fev std 14:5 5:914:9 3:615.1 4.1cpu std 2:2 0:82:9 0:82.8 0.825 Table 4: Results of analysis of the citrate signal with VARPRO.peak k fk std (kHz)dk std (kHz) ak std (a.u.)0k ( ) a1-0.0949 0.0004 -0.0222 0.0032 9.70 1.3002-0.1069 0.0003 -0.0251 0.0040 10.42 1.9703-0.1140 0.0005 -0.0111 0.0023 1.55 0.05 -98.744-0.1296 0.0001 -0.0126 0.0006 19.51 0.6435.605-0.1317 0.0001 -0.0129 0.0006 19.51 0.64 -35.606-0.1472 0.0005 -0.0111 0.0023 1.55 0.0598.74tted zero order phase ( ) 225.6297 25.8974tted begin time (ms) 2.4339 0.5528a 0k =0k 180= expresses the phase in degrees26 Table 5: Results of analysis of the citrate signal with the new algorithm.peak k fk std (kHz)dk std (kHz) ak std (a.u.)0k ( ) a1-0.0949 0.0004 -0.0227 0.0032 9.93 1.2702-0.1068 0.0003 -0.0240 0.0037 9.95 1.8303-0.1137 0.0001 -0.0117 0.0025 1.55 0.05 -98.744-0.1296 0.0001 -0.0127 0.0005 19.50 0.6335.605-0.1317 0.0001 -0.0127 0.0005 19.50 0.63 -35.606-0.1477 0.0001 -0.0117 0.0025 1.55 0.0598.74tted zero order phase ( ) 221.1242 25.2350tted begin time (ms) 2.3313 0.5380a 0k =0k 180= expresses the phase in degrees27 Captions for guresFigure 1: Frequency domain representation of simulated 31P MRS signal, the standard devia-tion of the noise is 15. See Table 1 for the used parameter values of the signal.Figure 2: In vivo spectrum of the prostate gland of a BPH patient. Peak 1 is choline, peak 2creatine and peaks 3 to 6 constitute the citrate multiplet.Figure 3: Graphical result of the analysis of the citrate signal with the new algorithm. Frombottom to top: the FT spectrum of the original signal, the individual Lorentzians (FT of ttedsinusoids) and the residual, which is the di erence between the original signal and the recon-structed signal.28 −200−1000100200300400500600700HzReal part of simulation signal 12345678